function fig = PlotRbSpinPolarization(~, polarizationData, field2freqCoeff)
    fittingRes = fittingData(polarizationData, field2freqCoeff);
    pwrList = linspace(0, max(polarizationData.pumpPower), 101);
    
    fig = figure('Name', sprintf('Rb Polarization Summary[%s]',datestr(datetime('now'))), 'WindowStyle', 'docked');
    ax1 = subplot(2, 2, 1);
    plot(ax1, polarizationData.pumpPower, polarizationData.width*field2freqCoeff/1e3, 'bd', ...
              pwrList, fittingRes.width.result(pwrList), 'r-' ); grid on; 
    xlim(ax1, minmax(polarizationData.pumpPower(:)')); xlabel(ax1, 'Pump Power (mW)'); ylabel(ax1, 'Peak Width (kHz)');
    legend(ax1, {'data', sprintf('f_0=%3.2fkHz\n k=%3.1f Hz/mW', fittingRes.width.f0, fittingRes.width.k*1e3)}, 'Location', 'best', 'Interpreter', 'tex');
%     legend(ax1, {'data', sprintf('f_0=%3.2fkHz\n k=%3.1f Hz/mW\n P_0=%3.1f mW', fittingRes.width.f0, fittingRes.width.k*1e3, fittingRes.width.p0)}, 'Location', 'best', 'Interpreter', 'tex');
    title(ax1, sprintf('\\Gamma_{20} =2\\pi\\times%3.2f kHz',  fittingRes.Gamma20_by2pi));
    
    ax2 = subplot(2, 2, 2);
    plot(ax2, polarizationData.pumpPower, polarizationData.hight, 'bd', ...
              pwrList, fittingRes.area.result(pwrList)./fittingRes.width.result(pwrList), 'r-'); grid on; 
    xlim(ax2, minmax(polarizationData.pumpPower(:)')); xlabel(ax2, 'Pump Power (mW)'); ylabel(ax2, 'Peak Hight (V)');
    title(ax2, sprintf('k_{R}=2\\pi\\times%3.2f Hz/mW', fittingRes.kR*1e3));
    
    ax3 = subplot(2, 2, 3);
    plot(ax3, polarizationData.pumpPower, polarizationData.width.*polarizationData.hight*field2freqCoeff/1e3, 'bd', ...
              pwrList, fittingRes.area.result(pwrList), 'r-'); grid on;
    xlim(ax3, minmax(polarizationData.pumpPower(:)')); xlabel(ax3, 'Pump Power (mW)'); ylabel(ax3, 'Area (V\cdot kHz)');
    legend(ax3, {'data', sprintf('P_0=%3.2f mW', fittingRes.area.P0)}, 'Location', 'best', 'Interpreter', 'tex');
    title(ax3, sprintf('\\Gamma_{10}=2\\pi\\times%3.2f kHz', fittingRes.Gamma10_by2pi));
    ax4 = subplot(2, 2, 4);
    shiftData = (polarizationData.bz0-polarizationData.bz0(1))*field2freqCoeff/1e3;
    plot(ax4, polarizationData.pumpPower, shiftData, 'bd', ...
              pwrList, fittingRes.shift.result(pwrList)); grid on;
    xlim(ax4, minmax(polarizationData.pumpPower(:)')); ylim(ax4, minmax(shiftData(:)')); 
    xlabel(ax4, 'Pump Power (mW)'); ylabel(ax4, 'Shift (kHz)');
    title(ax4, sprintf('Slope=%3.2f Hz/mW', fittingRes.shift.slope*1e3));
end

function fittingRes = fittingData(polarizationData, field2freqCoeff)
    %%
    [powerData, widthData] = prepareCurveData( polarizationData.pumpPower, polarizationData.width*field2freqCoeff/1e3 );
    ft = fittype( 'poly1' );
%     opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
%     opts.Display = 'Off';
%     opts.Lower      = [0 0 0 0];
%     opts.StartPoint = [widthData(end)/powerData(end) widthData(1) widthData(1)*mean(powerData) mean(powerData)];
%     opts.Upper      = [inf inf inf inf];
%     [width.result, width.gof] = fit( powerData, widthData, ft, opts);
    [width.result, width.gof] = fit( powerData, widthData, ft);
    width.f0 = width.result.p2;
%     width.p0 = width.result.q1;
    width.k = width.result.p1;
    fittingRes.width = width;
    
    fittingRes.Gamma20_by2pi = width.f0;
    fittingRes.kR = 2.0*width.k;
    %%
    [powerData, areaData] = prepareCurveData( polarizationData.pumpPower, polarizationData.width.*polarizationData.hight*field2freqCoeff/1e3 );
    
    % Set up fittype and options.
    ft = fittype( 'rat11' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Lower =      [0             0  0];
    opts.StartPoint = [max(areaData) 0  mean(powerData)];
    opts.Upper =      [Inf           0  max(powerData)];
    
    % Fit model to data.
    [area.result, area.gof] = fit( powerData, areaData, ft, opts );
    area.amplitude = area.result.p1;
    area.P0 = area.result.q1;
    fittingRes.area = area;

    fittingRes.Gamma10_by2pi = width.k*area.P0;
    
    %%
    [powerData, shiftData] = prepareCurveData( polarizationData.pumpPower, (polarizationData.bz0-polarizationData.bz0(1))*field2freqCoeff/1e3 );
    len = length(shiftData);
    ft = fittype('poly1');
    excludedPoints = excludedata( powerData, shiftData, 'Indices', floor(len/2):len );
    opts = fitoptions( 'Method', 'LinearLeastSquares' );
    opts.Lower = [-Inf 0];
    opts.Upper = [Inf 0];
    opts.Exclude = excludedPoints;
    % Fit model to data.
    [shift.result, shift.gof] = fit( powerData, shiftData, ft, opts );
    shift.slope = shift.result.p1;
    fittingRes.shift = shift;
end

